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Abstract. We develop a semi-classical method to simulate the motion of atoms in a dissipative optical 
lattice. Our method treats the internal states of the atom quantum mechanically, including all nonadiabatic 
couplings, while position and momentum are treated as classical variables. We test our method in the one- 
dimensional case. Excellent agreement with fully quantum mechanical simulations is found. Our results are 
much more accurate than those of earlier semi-classical methods based on the adiabatic approximation. 

PACS. 32.80.Pj Optical cooling of atoms; trapping - 03.65.Sq Semiclassical theories and applications 



1 Introduction 



One of the most spectacular achievements in the field of 
laser cooling is the discovery of cooling below the Doppler 
limit in optical lattices, so called Sisyphus cooling 1 . An 
optical lattice is a standing wave of laser light, forming 
a periodic light-shift potential for atoms moving in the 
laser field 12113 • In the optical lattices used for cooling 
the frequency of the lasers are tuned close to an atomic 
resonance. The atoms thus undergo cycles of absorption 
followed by spontaneous emission. Under the right exper- 
imental conditions, the spontaneous emission causes an 
overall loss of kinetic energy of the atoms, i.e., cooling. 

Optical lattices are also widely used in Bose-Einstein 
condensation experiments A and for quantum state ma- 
nipulation la]. These lattices are tuned far from atomic 
resonances, in order to avoid spontaneous emission which 
would destroy the coherence of the condensate. Therefore 
these far detuned lattices do not provide any cooling. 

The name Sisyphus cooling comes from the first the- 
oretical model for the process [S||7]. This model is based 
on optical pumping between the magnetic sublevels of the 
light shifted atomic ground state. However, at least in its 
original form it relies on a number of simplifying assump- 
tions, such as a semi-classical approximation, spatial av- 
eraging, and a simplified level structure (a ground state 
with angular momentum J g = 1/2, and an excited state 
with angular momentum J c = 3/2). Whereas this model 
correctly predicts some qualitative features of cooling in 
optical lattices, it is too crude to provide an overall quan- 
titative agreement. Instead, a number of more advanced 
theoretical methods have been developed. The most accu- 
rate of these is the Monte-Carlo wavefunction technique 
|S] , a fully quantum mechanical method based on stochas- 
tic wavefunctions. 



In this paper, we develop and test a new semi-classical 
method for simulating the motion of atoms in a near- 
resonant optical lattice. The most important approxima- 
tion of our method is that the position and momentum of 
the atoms are treated as classical variables. Other approx- 
imations include a classical treatment of the light field, 
and adiabatic elimination of excited states of the atoms, 
but otherwise we make as few approximations as possible. 
In particular, the internal states are treated quantum me- 
chanically, allowing for any kind of coherent superposition 
between magnetic sublevels. 

Even though more exact fully quantum mechanical 
theoretical methods exist, semi-classical methods are valu- 
able, partly because they are less demanding numerically, 
but also because they provide a simpler conceptual frame- 
work in which it is easier to formulate an intuitive picture 
of e.g. the mechanisms involved in the cooling process. Up 
to now, all semi-classical methods for laser cooling in op- 
tical lattices have been based on atoms that are pumped 
between definite internal states as they move through the 
lattice. To this end a basis of so-called adiabatic states, di- 
agonalizing the light-shift potential at every position, has 
been used instead of the diabolic basis of the magnetic 
substates UJ. Coherences between adiabatic states have 
not been included in the description, and neither have 
so-called nonadiabatic couplings arising from the position 
dependence of the adiabatic basis. Thus, the motion of 
the atoms is described by purely classical equations, al- 
beit the various potentials, pumping rates and diffusion 
coefficients have been derived from a quantum-mechanical 
origin. These adiabatic semi-classical methods reproduce 
some of the qualitative features of Sisyphus cooling, e.g. a 
linear relation between temperature and irradiance at high 
irradiances |l()j . However, we show that even at very high 
irradiances the slope of this linear dependence does not 
agree with fully quantum-mechanical simulations. At the 
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lower irradiances relevant to most experiments the adia- 
batic semi-classical method deviates even more severely 
from the fully quantum- mechanical results. Both these 
problems are solved by the nonadiabatic semi-classical ap- 
proach. 



2 Theory 

In this section we develop the basic semi-classical equa- 
tions of motion, on which our simulations are based. For 
generality the theory is developed in three dimensions. 
The angular momenta of the ground and excited states 
of the lattice transition are denoted by J g and J e respec- 
tively, and the corresponding magnetic quantum numbers 
are M g and M c . Although the derivation is more general, 
we shall in the end apply the theory to the case J c = J g +1- 
Also, the light field £(r) creating the lattice could take 
different forms, but will in the end be assumed to have a 
lin_Llin configuration in one, two or three dimensions 3 . 
That is, the lattice is created by the interference pattern of 
light fields, forming lattice sites with alternating a + and 
a~ polarizations. 

We start from the optical Bloch equations for an atom 
in a standing wave laser field They can be derived 
under very general conditions, and represent for practical 
purposes an exact fully quantum mechanical description 
of atomic motion in an optical lattice. Our first important 
approximation is that the population of the excited state 
is sufficiently low to allow its adiabatic elimination. The 
condition for this is that the saturation parameter 



so = 



f2 2 /2 



A 2 + T 2 /A 



< 1. 



(1) 



Here A is the detuning from resonance, F the natural 
width of the excited state, and Q is the Rabi frequency 1 . 
The details of the adiabatic elimination of the excited 
states can be found e.g. in Ref. (TJ. This approximation 
is an important simplification, since it reduces the master 
equation for the full density matrix, to an equation for the 
(2 J g + 1) x (2 J g + 1) density matrix a of the ground states. 
The resulting equation for the evolution of a reads 



H,<7 



ih a 



sp ■ 



(2) 



The first term on the right-hand side of this equation rep- 
resents the Hamiltonian part of the evolution. The second 
term represents the non-Hermitian evolution due to spon- 
taneous emission. The Hamiltonian contains the kinetic 
term and the light-shift potential, 



H=^ 

2m 



HA'A{r), 



(3) 



We use the Rabi frequency based on the total laser field. 
This is the same convention as was used, e.g., in Ref. 0. Some- 
times the Rabi frequency is instead on the laser irradiance per 
beam, which for a one-dimensional lin_Llin configuration is half 
the total irradiance. 



where p is the momentum operator of the atom, r its 
position, A' = Asq/2, and the operator A(r) is given by 



A(r) = 



■e(r) 



(4) 



Here d + is an operator that promotes an atom from the 
ground to the excited state, while d~ = (d + )^ is respon- 
sible for the reverse process. In the basis of circular polar- 
ization vectors 



£±i=T^(x±y), £ = z, 



(5) 



they have simple expressions in terms of Clebsch-Gordan 
coefficients 



d+ = (J c Mc | JglM g q) = [d 



(6) 



In the basis of the magnetic substates M g the operator 
A(r) is represented a matrix A(r). 

For the simple model atom with J g = 1/2 and J c = 3/2 
A(r) is a diagonal matrix. However, most atoms of interest 
have a more complicated level structure, including non- 
diagonal couplings in the potential. Therefore previous 
semi-classical methods have used an adiabatic basis, where 
the atomic states are the eigenstates of A{r). Whereas 
A(r) is diagonal in the adiabatic basis, the position de- 
pendence of the basis gives rise to nonadiabatic couplings 
between adiabatic states. In the adiabatic approximation 
these couplings are neglected. In our method we keep all 
off-diagonal couplings. The results arc then independent 
of the basis used, and the simplest choice is to stay with 
the magnetic levels Af g , the diabolic basis. Since this basis 
is the same for all r, all couplings are included in A(r), 
and their functional form can be calculated analytically 
for a given laser configuration. 

The second term on the right-hand side of Eq. (0) 
contains processes associated with spontaneous emission. 
Writing the matrix elements of a in the position represen- 
tation, (r\a\r')= er(r,r'), its form is 

r' 3-T' 

°(r,r% p =- — [A{r)<j{ry) + a{ry)A(r')] + — 

f <W k J2 Bl(r)c- ik r a(r, r')e ik r ' B e (r'), 

(7) 



e±k 



where r' = _Tsq/2. The matrices B e (r) are given by 



Be(r) 



(8) 



Hence, B\ represents the excitation of an atom by the 
laser field, and its subsequent return to the ground state 
via spontaneous emission of a photon with polarization e. 
The factors exp(ifc • r) represent the atomic recoil from 
a spontaneously emitted photon with wave vector k. The 
integration is over the directions of the emitted photon, 
and the summation is over any basis spanning the allowed 



S. Jonsell et al.: A nonadiabatic semi-classical method for dynamics of atoms in optical lattices 



3 



polarization vectors. The recoil momentum of the atomic 
transition is j»r = hk^ = h\k\. 

Our goal is to approximate Eq. (J2J) by a semi-classical 
equation where every atom has a definite position and 
momentum, i.e. every atom follows a trajectory in phase 
space. This is of course not allowed in quantum mechanics, 
because of the uncertainty principle. Hence, a quantum 
mechanical phase space cannot be defined, but it is still 
possible to introduce a "coarse grained" version of phase 
space through the Wigner function 



W(r,p,t) = ^/ dM ( r + | 



r _|^ c -ip-Wi. (9) 



In this work the Wigner function is a matrix with dimen- 
sion 2J g + 1. The Wigner transformation of Eq. J2J be- 
comes 



d_ 

dt 



p_ 

rn 



V r )W(r,p,t) 
A' 



fi3 

r 



dqe^ r ^[W(r,p+^,t)A(q) 
-A(q)W(r,p-^t)] 
dqj**/*[W(r,p+Z,t)A(q) 
+ A(q)W(r,p-^,t)] 



x Bl(q')W(r,p + hk 



Here A and B e are the Fourier transforms 
A(q) = J dre-' lq r / h A(r), 



,t)B e (q). 
(10) 

(11) 
(12) 



No approximation has been made in going from Eq. 
(J5J to Eq. (|10fl . the Wigner transformation is just another 
representation of the same physics. Now, we introduce the 
semi-classical approximation. According to this approxi- 
mation the momentum distribution varies smoothly and 
slowly over typical momentum transfers q in Eq. 
Since A(r) and B e (r) have the same periodicity as the 
laser field, i.e. A = 27r/fcrt, Eqs. (|ll(l and Ijl2|l show that the 
typical size for q is the recoil momentum. Thus, the semi- 
classical approximation assumes that the momentum dis- 
tribution changes little for emission/absorption of a single 
photon. As long as the atomic momenta are several recoil 
units large, and the effects of quantization of the atomic 
states are small, this approximation can be expected to 
work well. 

Invoking the semi-classical approximation we can make 
a second-order Taylor expansion around p of the Wigner 



distribution 

W{r,p + q, t) ~W{r,p, t) + q ■ V p TU(r,p, t) 

+ \{<1- V p ) 2 W(r,p,t). (13) 

Using this expansion it is possible to replace A and B e by 
their counterparts in position space. The resulting equa- 
tion for the semi-classical Wigner function, which now can 
be interpreted as a phase-space distribution, is 



i=l 



r' 



{W(r,p,t),A(r)} 



r' J2 Bl(r)W(r,p,t)B q (r) 

q=0,±l 
3 



HA' 



^{d Pi W(r,p,t),diA(r)} 

3 

Y,[9piW{r,p,t),diA(r)] 



.HP 



i=l 

3 



i— J2 J2^ B l^ w ( r ^^ B ^ 



q=0,±l i=l 



Bj(r)0 w W(r,p,t)$B,(r) 



h 2 A 



i=i j=i 

h 2 r' 3 3 



1G 

h 2 r 



i=l 3 = 1 
2rt 33 
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q=0,±l i=l j=l 

- 2diBl (r)d Pi d P] W(r, p, t)d B q (r) 
+ Bl(r)d Pi d Pj W(r,p,t)d i d j B q (r)] 



3 



<J=0,±1 i=l 



(14) 



In this equation we use the short-hand notation di = 
d/dri, d Pi = d/dpi, where i — x,y,z are the Cartesian 
coordinates. The constants T]i, q come from the integration 
over the direction of the spontaneously emitted photon, 
and are given by 

Vx,±i = Vy,±i = 3/4, r) zfi = 1/2 

Vx,0 = Vy,o = 1, Vz,±l = 1- (15) 

Although the equation is somewhat lengthy, it is possible 
to give physical interpretations to its terms. The left hand 
side is simply the kinetic term, i.e. the full derivative d/dt. 
On the right, the terms where W(r,p, t) appear without 
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any derivative represent transfer of population between 
states, either by couplings from non-diagonal terms of the 
light-shift potential A'A(r), or by optical pumping. The 
terms containing d Pi W(r,p, t) describe the motion of the 
atoms due to forces from light-shift potential and the ra- 
diation pressure. Terms containing second derivatives of 
both W(r,p,t) and first or second derivatives of A{r) 
or B e (r) describe the momentum diffusion due to fluc- 
tuations in the number of photons absorbed. Finally, the 
term containing dp.W(r,p,t), but no other derivatives, 
contains the momentum diffusion due to the recoil kick 
from spontaneously emitted photons. 

Equation (|14(l is the most complete semi-classical ap- 
proximation for the time-dependent distributions of atoms 
in r and p space. It is classical in the sense that the 
atoms are assumed to be particles with definite positions 
and momenta. The internal states, however, are treated 
fully quantum mechanically, including all off-diagonal cou- 
plings and coherences. It is thus not possible to assign an 
atom to a definite internal state, nor is it described as a 
classical probability distribution over the different inter- 
nal states, but as a quantum-mechanical superposition of 
internal states. 

In order to solve Eq. (|14|l we recast it into a Langevin- 
type equation. That is, instead of calculating distributions 
of atoms, we shall calculate phase-space trajectories x(t) 
and p{t) of individual atoms. In doing this, we still want to 
keep the quantum mechanical description of the internal 
states. That is, the probability distribution of an atom is 

W{r, p, t) = w(t)6 (r - r(t)) 8 (p - p(t)) . (16) 

Here w(t) is a matrix of dimension 2J g + 1 containing 
the internal-state density matrix of the atom at time t. 
Inserting this form into Eq. I|14(l , and integrating over po- 
sition and momentum, the evolution equation for w(t) is 
obtained 

w(t) =iA'[w(t),A(r)]-^{w(t),A(r)} 

+ T' B\{r)w{t)B q {r). (17) 

g=0, ±1 

Here and below, we use the simplified notation r for r(t) 
and p for pit). It is, however, important to understand 
that these are now time-dependent functions represent- 
ing position and momentum of a single atom, which are 
conceptually very different from the variables in Eq. (|14|) . 
Using that (x) = Tr{a;u;} etc., we derive the equations 
for the evolution of x and p (see, e.g., \V2\) 

* = (18) 

TO 

V = f(t)+x(t). (19) 

Here f(t) is a force and x(t) is a fluctuating force with 
the properties 

(*(*)) = 0, (Xi(t)Xj(t')) = 2 A, (t)5(t - t'). (20) 



The force is given by 

f i (t)=-HA'Tr{d i Air)w(t)} 

-iy E Tr{[B q ir)d t Bl(r) 

9=0, ±1 

-ftfl,(r)Bt(r)]«,(i)}. (21) 

The first term above is the force arising from the second- 
order light-shift potential, while the second term is the 
radiation pressure. The diffusion coefficient is 

!>«(*) =««^^ E Vi, q Tr{B q (r)Bl(r)w(t)} 
?=o,±i 

+ ^-) E H[Wr)d 3 Bi(r) 

x J ' q— 0,±1 

+ d 3 B q {r)d t B\{r)\w{t)}. (22) 

The first term arises from the recoil from photons spon- 
taneously emitted in random directions, while the second 
term is connected to fluctuations in the radiation pressure. 
The latter term is in general anisotropic. 



3 Numerical implementation 

We simulate the equations i|17|) . i|18[) and (|19(l in one di- 
mension. The laser field has the form 

£(z) = cos(fcRz)e_i — isin(fcR,z)£ + i. (23) 

At the start of every time step the system is in a pure 
quantum mechanical state. For every time step w,z, and 
p are evolved using a second-order Runge-Kutta method. 
The fluctuating force x(i) is included as a term 

rV2D dt, (24) 

where r is a random number with zero average and unit 
variance. This term only needs to be evaluated once every 
time step [T3] . 

At the end of a time step, the system will not be in 
a pure state anymore. Its internal-state density matrix w 
can, however, be decomposed into 2 J g + 1 pure states 

2J g + l 

w= E ( 25 ) 

i=l 

The coefficient are the eigenvalues, and \<Pi) the cor- 
responding eigenvectors, of w. Since w is a density ma- 
trix, the eigenvalues satisfy the properties > and 

J2i=\ +1 — 1) an d can be interpreted as classical proba- 
bilities of the different states ^31- Hence, at the end of 
each time step the system has the probability A.; to make 
a "jump" into the pure state |<2>j). Even though a density 
matrix in general has an infinite number of decomposi- 
tions into pure states, the decomposition above is unique 
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in the sense that it is the only one into a set of linearly 
independent pure states. 

For numerical efficiency the eigenvalues were obtained 
using first-order perturbation theory, which is sufficiently 
exact if dt is short enough. In practice, one of the eigenval- 
ues will be very close to one, while the others are small or 
zero. Thus, one can interpret the system as either staying 
in the same state, or jumping to a new state. When the 
eigenvalues obtained by perturbation theory indicate that 
the system makes a jump, the accuracy is increased by a 
full diagonalization of w. The expense in computer time 
for this improvement is modest, since jumps are compar- 
atively rare. 



4 Results 

In our simulations we used the parameters for the D2 line 
in cesium, i.e. J g — 4, J c = 5, and natural width r/2ir = 
5.2227 MHz, and recoil energy E R = 1.3692 x 1CT 30 J 
|15| . The diagonal elements of the diabatic potential for 
this transition are displayed in Figure ^ We first investi- 
gated the steady-state momentum distributions. For po- 
tential depths H\A'\ > 200-Er the samples contained 5000 
atoms, and were iterated for the time 2500/ J". To im- 
prove statistics the momentum distribution was averaged 
over the last 1000/J 1 ' of the evolution time. For low poten- 
tial depths convergence is slower. Therefore we used 20000 
atoms for HA' < 200_Er, and the evolution time 5000/7"', 
with averaging over the last 2000/-T'. For all runs the time 
step was dt = 0.025/ T", and the initial state a spatially 
uniform distribution with temperature of 10 /iK. 



results are compared to a full-quantum simulation using 
the Monte-Carlo wave function method [B] • The two meth- 
ods are in excellent agreement. The relative difference is 
at most about 20%. It is not clear how much of this de- 
viation can be attributed to the fundamental difference 
between the two methods, and how much is due to e.g. sta- 
tistical uncertainties or other numerical errors. For deep 
potentials both methods give the same linear slope, al- 
though with a slight offset. The agreement continues all 
the way down through decrochage, i.e. the point where the 
curve turns around and starts to increase again for small 
potential depths, although statistical fluctuations in the 
full-quantum data make comparisons more difficult here. 

It also evident from Figure |21 that the present method 
is a substantial improvement of the adiabatic method used 
in Ref. [Hj- The methods do not even agree at large po- 
tential depths, where one would expect the nonadiabatic 
corrections to become small. Improving upon this method 
by including non-diagonal diffusion terms (for details see 
Ref. 0) does not substantially change the situation. We 
note that even in the limit of vanishing nonadiabatic cor- 
rections our method differs from that in Ref. [H] by al- 
lowing for coherences between the internal states. In the 
adiabatic basis the potential does not induce any coher- 
ences between internal states, but such coherences are still 
induced by optical pumping. 
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Fig. 1. Diagonal elements of the diabatic potential for the 
J g = 4 — » J = 5 transition. Each curve corresponds to a mag- 
netic sublevel M g of the ground state. Curves corresponding 
to ± | A- J g share the same color coding, and differ only by the 
phase 7r/2. States with M s even (odd) are represented by solid 
(dashed) curves. 



Results for (p 2 ) as a function of potential depth HA' / E R 
for a detuning A — — 10J" are displayed in Figure Our 
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Fig. 2. Semi-classical results for (p 2 ) (circles) compared to full 
quantum results (squares). For comparison we also show results 
based on the adiabatic approximation calculated similarly to 
the method used in Ref. (crosses), and the same method 
improved by including also non-diagonal diffusion coefficients 
(triangles). The detuning is A = — 10J\ 



The semi-classical method also makes it possible to fol- 
low the motion of a single atom as it moves through the 
lattice. In Figures [3] and Q] we show the position, momen- 
tum, energy and internal state distribution as a function 
of time for a single atom in optical lattices with detun- 
ings A — — 10r, and potential depths K\A'\ — 150£r and 
H\A'\ = IOOOE'r respectively. The energy was calculated 
as the sum of the kinetic energy and light-shift potential, 
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dynamics of atoms in optical lattices 




rt 



Fig. 3. Position, momentum, energy, and internal state pop- 
ulations as a function of time for a single atom moving in an 
optical lattice. The potential depth is fr\A'\ = 150i?R, and the 
detuning A = — KIT. The internal states have the color coding 
from Figure Q 



i.e., 

E=£- + hA'Tr{A(z)w(t)}. (26) 
2m 

The ratio between the potential, pumping and diffu- 
sion terms in Eq. (|14|) depends on A/T only, and is hence 
the same in both graphs. The only difference lies in the in- 
ertial term p/mdi. Increasing \A'\, while keeping the ratio 
A/r constant, is equivalent to increasing the mass m by 
the same factor. This can be seen comparing the graphs, 
since the atom is less mobile in Figure 0J 

At both potential depths the atom shows, after an ini- 
tial cooling phase, a high degree of localization. While lo- 
calized the atom populates mostly the extreme magnetic 
states M g = ±J g . The energy is more or less constant, 
fluctuating around half the potential depth. The ampli- 
tudes of the oscillations in momentum and position vary 
somewhat due to diffusion, but tend to stay within certain 
bounds as long as the atom remains in the same poten- 
tial well. We cannot see any clear trend towards smaller 




rt 

Fig. 4. Same as Figure [3] but for a deeper potential, h\A'\ — 
IOOOSr . 



oscillation amplitudes while the atom remains trapped in 
a site, i.e., we see no local cooling. 

The periods of localization are interrupted by brief 
phases where the atom acquires enough energy to travel 
over many potential wells, before once again getting lo- 
calized. These excursions are most prominent at lower po- 
tential depths. The periods when the atom is untrapped 
are associated with abrupt changes of the internal state of 
the atom, usually from odd to even magnetic states. (The 
light-shift potential only induces odd-odd and even-even 
couplings between magnetic states. Thus any pure quan- 
tum mechanical state is a superposition of only odd or only 
even magnetic states.) During all periods of localization 
the atom is in a state with similar internal-state distribu- 
tion and energy. Even when the energy sometimes drops 
below this stationary value the atom is soon returned to 
the same state. 

These results are in qualitative agreement with our 
earlier conclusion that Sisyphus cooling, especially at low 
potential depths, works through a transfer of atoms be- 
tween a hot and a cold mode The cold mode has 
a momentum distribution, with a width that does not 
change over time. This mode corresponds to the popu- 
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lation of atoms in the trapped state. The cooling process 
is in effect a transfer of atoms from the untrapped to the 
trapped state. 

In Fig. we compare the semi-classical approxima- 
tion to the time evolution of the momentum distribu- 
tion D(p) — AN{p)/dp (where N(p) is the number of 
atoms with momentum p) to the results in for \A'\ = 
130-Er. The bimodality of the distribution is very clear 
also in the semi-classical results, and the agreement with 
the quantum-mechanical results is very good. The distri- 
bution of the hot mode is identical to within statistical 
uncertainties. This shows that the physics of untrapped 
atoms, including their rate of transfer to trapped states, 
is well described by our semi-classical method. The semi- 
classical method gives a slightly more narrow cold mode, 
in agreement with the results in Fig. [21 
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Fig. 5. Time evolution of the momentum distribution for a 
potential depth \A \ — 130-Er. The starting temperature was 
50/iK. The black curve shows the semi-classical results, while 
the red curve shows results of a fully quantum mechanical sim- 
ulation. 



5 Discussion 

We have developed a semi-classical method to simulate the 
dynamics of atoms in optical lattices. Our results for the 
average momentum distribution of the atoms, including 
its time dependence, agree excellently with those of the 
fully quantum mechanical method. To achieve an accurate 
description it is necessary to include both populations of 
and coherences between the internal states of the atom. 
The external degrees of freedom may, at least in some 
situations, be described classically, i.e., as particles with 
definite positions and momenta. 

The semi-classical approximation was introduced as a 
second order Taylor expansion in p/pn of the Wigner func- 
tion. According to our results (p) Tms > ipn, and hence this 
expansion should be a fairly good approximation. Never- 
theless, there are some situations where the semi-classical 
description must necessarily break down. One is when ef- 
fects from the quantization of bound states are important. 



Such effects will be most prominent when the atoms are 
localized near the bottom of the potential wells. Another 
is the coherent splitting of a wave packet. If the atomic 
wavefunction is, e.g., partially transmitted to the next po- 
tential well, the semi-classical method will describe this as 
a classical probability (some atoms are transmitted, some 
are not), while any coherence effects between the two parts 
of the wave packet will be lost. 

The conceptual simplicity of the semi-classical descrip- 
tions makes it a useful aid in visualizing complex physical 
processes. It is also a flexible tool, which is relatively easy 
to adapt to different physical situations. In the near future 
we plan to extend the method to double optical lattices 
[T7|. Further studies of the cooling process, e.g. to deepen 
the understanding of the bimodal velocity distributions 
observed in experiment and full quantum simulations, are 
underway. 
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